
# Loading -----------------------------------------------------------------

setwd("C:/Users/elies/Dropbox/Climate fairness/06 Empirics/05 DCP data/07 Replication NComms")

library(haven)

deppack <- c("glmnet", "lars", "Matrix", "quadprog", "glinternet", "igraph", "sandwich",
             "lmtest", "stats", "graphics", "utils", "limSolve", "arm")
for(pkg in deppack){
    eval(bquote(library(.(pkg))))
}

source("01 Code/05 amie support functions.R")
## this file contains the entire FindIt package
## with slight adjustments made to two functions

dt <- read_dta("02 Data/pooled_scenlevel_touse.dta")
dt <- subset(dt, !is.na(dt$own_cstpth_cj))

# Preparing data ----------------------------------------------------------

dt$own_hhcsts_l <- factor(dt$own_hhcsts_cj,
                          levels = 1:4, ordered = T,
                          labels = c("low", "medium", "high", "vhigh"))
dt$other_hhcsts_l <- factor(dt$other_hhcsts_cj,
                            levels = 1:4, ordered = T,
                            labels = c("low", "medium", "high", "vhigh"))
dt$own_mitig_l <- factor(dt$own_mitig_cj_5levels,
                         levels = 1:5, ordered = T,
                         labels = c("0-19", "20-39", "40-59", "60-79", "80-100"))
dt$other_mitig_l <- factor(dt$other_mitig_cj_5levels,
                           levels = 1:5, ordered = T,
                           labels = c("0-19", "20-39", "40-59", "60-79", "80-100"))
dt$own_cstpth_l <- factor(dt$own_cstpth_cj,
                          levels = c(2,1,3), ordered = T,
                          labels = c("increasing", "constant", "decreasing"))
dt$other_cstpth_l <- factor(dt$other_cstpth_cj,
                            levels = c(2,1,3), ordered = T,
                            labels = c("increasing", "constant", "decreasing"))

# Analysis ----------------------------------------------------------------

db <- "#3e647d"
lb <- "#008bbc"
refcat <- rep(4, 21)

mod.a <- CausalANOVA(formula = choice_cj ~ own_hhcsts_l + other_hhcsts_l +
                         own_mitig_l + other_mitig_l + own_cstpth_l + other_cstpth_l,
                     data = dt, nway = 2, cluster = dt$caseid)

ylabsA <- c("Own Low x Other Low", "Own Medium X Other Low", "Own High X Other Low",
            "Own V High X Other Low", "Own Low X Other Medium", "Own Medium X Other Medium", 
            "Own High X Other Medium", "Own V High X Other Medium", "Own Low X Other High",
            "Own Medium X Other High", "Own High X Other High", "Own V High X Other High",
            "Own Low X Other V High", "Own Medium X Other V High", "Own High X Other V High",
            "Own V High X Other V High")
colsA <- c(lb, rep(db, 4), lb, rep(db, 4), lb, rep(db, 4), lb)

png(filename = "03 Results/Fig5.png", width = 10, height = 8, units = "in", res = 400)
plot(mod.a, type = "AMIE", c("own_hhcsts_l", "other_hhcsts_l"), space = 15,
     maintitle = "Average Marginal Interaction Effects",
     ylabs = rev(ylabsA), col = colsA)
dev.off()
